#pragma once

#define _USE_MATH_DEFINES
#include <cmath>
#include <algorithm>

#include <QColor>
#include <QSet>
#include <QVector>
#include "float.h"

#include "geometry/Vector.h"

inline unsigned int fact(unsigned int n){
    unsigned int i, p=1;

    for(i = 2; i <= n; i++)
        p *= i;

    return p;
}

// Utility Macros
#define Max(a,b) (((a) > (b)) ? (a) : (b))
#define Min(a,b) (((a) < (b)) ? (a) : (b))
#define MaxOf(a,b,c) (Max(a, Max(b,c)))
#define Mod(x,m) ((x % m + m) % m)
#define RANGED(min, v, max) ( Max(min, Min(v, max)) ) 
#define RANGE(i, min, max) (  ((i >= min) && (i <= max)) ? 1 : 0)
#define SIGN(i) ((i >= 0) ? (1) : (-1))
#define SET_HAS(SET, OBJECT) ((SET.find(OBJECT) != SET.end()) ? 1 : 0)
#define MAP_HAS(MAP, OBJECT) ((MAP.find(OBJECT) != MAP.end()) ? 1 : 0)
#define PREV(i, N) ((i + N-1) % N)
#define NEXT(i, N) ((i + 1) % N)
#define SWAP(x, y, T) do { T temp##x##y = x; x = y; y = temp##x##y; } while (0)
#define AROUND(x, target, threshold) ( (abs(x) - abs(target) < threshold) ? 1 : 0)
#define RADIANS(deg)    ((deg)/180.0 * M_PI)
#define DEGREES(rad)    ((rad)/M_PI * 180.0)


// Array operations
#include <algorithm>
#define MaxElement(v) (*max_element(v.begin(), v.end()))
#define MinElement(v) (*min_element(v.begin(), v.end()))
#define DivideVector(a,value) for(int i = 0; i < a.size(); ++i) a[i] /= value;
#define Sum(v)  (std::accumulate(v.begin(), v.end(), 0.0))
#define Avg(v)  (Sum(v) / v.size())

// Rodrigues' rotation
#define ROTATE_VEC(v, theta, axis) (v = v * cos(theta) + cross(axis, v) * sin(theta) + axis * dot(axis, v) * (1 - cos(theta)))
#define ROTATED_VEC(v, theta, axis) (v * cos(theta) + cross(axis, v) * sin(theta) + axis * dot(axis, v) * (1 - cos(theta)))

// transformation
inline Vec3d reflect(Vec3d v, Vec4d plane)
{
	Vec3d n(plane[0], plane[1], plane[2]);	// normal
	Vec3d pp = -plane[3] * n;				// point on plane
	double dis = dot( pp - v, n);			// distance to plane

	return v + 2 * dis * n;
}

inline bool onPositiveSide(Vec3d v, Vec4d plane)
{
	Vec3d n(plane[0], plane[1], plane[2]);	// normal
	Vec3d pp = -plane[3] * n;				// point on plane
	double dis = dot( pp - v, n);			// distance to plane

	return dis < 0;
}

// Math
inline double round(double r) {
	return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}

// Random
inline double uniform_double(double a = 0.0, double b = 1.0)
{
	double len = b - a;
	return ((double)rand()/RAND_MAX) * len + a;
}

inline int uniform_int(int a, int b)
{
	double r = uniform_double(a, b);
	return int(round(r));
}

inline QSet<int> subset_indices(int m, int n)
{
	QSet<int> subset_ids;
	// m < n
	if (m > n-1) return subset_ids;

	int mm = (m <= n/2) ? m : n - 1 - m;
	while(subset_ids.size() < mm)
		subset_ids.insert(uniform_int(0, n-1));

	// The complement
	if (m > n/2)
	{
		QSet<int> all_ids;
		for (int i = 0; i < n; i++) all_ids.insert(i);

		return all_ids - subset_ids;
	}
	else
		return subset_ids;
}

inline void RandomBarycentric(double * interp)
{
	interp[1] = uniform_double();
	interp[2] = uniform_double();

	if(interp[1] + interp[2] > 1.0)
	{
		interp[1] = 1.0 - interp[1];
		interp[2] = 1.0 - interp[2];
	}

	interp[0] = 1.0 - (interp[1] + interp[2]);
}


// Coordinates transfer
inline Vec3d sphericalToEuclidean(double r, double theta, double phi)
{
	return r * Vec3d(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
}
